#include <iostream>
#include <cmath>
#include <climits>
#include <cstdlib>
#include <ctime>
#include <cstdio>

/* Ovaj define ukljucuje ispis koordinata svih uzoraka */
//#define MAMC_SAMPLE_LOGGING

using namespace std;

typedef long double LD;
typedef unsigned long long ULL;

ULL sampling = 0;
int _miserSteps = 2;

enum IntegrationMethod
{
	IM_Uniform,
	IM_MonteCarlo,
	IM_AdaptiveMC
};

typedef struct
{
	LD (*mathFunction1D)(LD);
	LD (*mathFunction2D)(LD, LD);
	LD lowerLimitX;
	LD upperLimitX;
	LD lowerLimitY;
	LD upperLimitY;
	LD multiplier;
	LD result;
} MAMC_TestCase;

/* Pokazivaci na trenutnu 1D i 2D funkciju */
LD (*_currentFunc)(LD) = NULL;
LD (*_currentFunc2D)(LD, LD) = NULL;

/* Obicna Monte Carlo integracija na koju se oslanjaju i ostale varijacije */
LD mcIntegrate(LD x1, LD x2, ULL samples, LD &mean, LD &ssum, LD &variance)
{
	LD sum = 0, var = 0;
	LD diff = x2 - x1;
	LD multi = diff / RAND_MAX;
	LD error;

	for(ULL i = 0; i < samples; i++)
	{
#ifdef MAMC_SAMPLE_LOGGING
		LD X = rand() * multi + x1;
		LD y = _currentFunc(X);
		printf("%.10Lf %.10Lf\n", X, y);
#else
		LD y = _currentFunc(rand() * multi + x1);
#endif
		sampling++;
		sum += y;
		var += y * y;
	}

	ssum = var;
	mean = sum / samples;

	error = (x2 - x1) * sqrt(((var / samples) - (sum / samples) * (sum / samples)) / samples);

	if(samples < 2)
		variance = 0.0;
	else
		variance = var / (samples - 1) - mean * mean;

	mean *= diff;
	return error;
}

/* Obicna Monte Carlo integracija u dve dimenzije */
LD mcIntegrate2D(LD x1, LD x2, LD y1, LD y2, ULL samples, LD &mean, LD &ssum, LD &variance)
{
	LD sum = 0, var = 0;
	LD diffX = x2 - x1;
	LD diffY = x2 - x1;
	LD multiX = diffX / RAND_MAX;
	LD multiY = diffY / RAND_MAX;
	LD error;

	for(ULL i = 0; i < samples; i++)
	{
#ifdef MAMC_SAMPLE_LOGGING
		LD X = rand() * multiX + x1, Y = rand() * multiY + y1;
		LD z = _currentFunc2D(X, Y);
		printf("%.10Lf %.10Lf %.10Lf\n", X, Y, z);
#else
		LD z = _currentFunc2D(rand() * multiX + x1, rand() * multiY + y1);
#endif

		sampling++;
		sum += z;
		var += z * z;
	}

	ssum = var;
	mean = sum / samples;

	error = (x2 - x1) * sqrt(((var / samples) - (sum / samples) * (sum / samples)) / samples);

	if(samples < 2)
		variance = 0.0;
	else
		variance = var / (samples - 1) - mean * mean;

	mean *= diffX * diffY;
	return error;
}

/* Uniformna integracija za jednu dimenziju */
LD uniformIntegrate(LD x1, LD x2, ULL samples, LD &mean, LD &ssum, LD &variance)
{
	LD sum = 0, var = 0;
	LD diff = x2 - x1;
	LD ldSamp = (LD)samples;
	LD error;

	for(ULL i = 0; i < samples; i++)
	{
#ifdef MAMC_SAMPLE_LOGGING
		LD X = i / ldSamp * diff;
		LD y = _currentFunc(X);
		printf("%.10Lf %.10Lf\n", X, y);
#else
		LD y = _currentFunc(i / ldSamp * diff);
#endif
		sampling++;
		sum += y;
		var += y * y;
	}

	ssum = var;
	mean = sum / samples;

	error = (x2 - x1) * sqrt(((var / samples) - mean * mean) / samples);

	if(samples < 2)
		variance = 0.0;
	else
		variance = var / (samples - 1) - mean * mean;

	mean *= diff;
	return error;
}

/* Uniformna integracija za dve dimenzije */
LD uniformIntegrate2D(LD x1, LD x2, LD y1, LD y2, ULL samples, LD &mean, LD &ssum, LD &variance)
{
	LD sum = 0, var = 0;
	LD diffX = x2 - x1;
	LD diffY = y2 - y1;
	LD ldSamp = (LD)samples;
	LD error;

	ULL dimSamp = sqrt(samples);
	samples = dimSamp * dimSamp;
	ldSamp = (LD)samples;
	LD ldDimSamp = (LD)dimSamp;

	for(ULL y = 0; y < dimSamp; y++)
	{
		for(ULL x = 0; x < dimSamp; x++)
		{
#ifdef MAMC_SAMPLE_LOGGING
			LD X = x / ldDimSamp * diffX + x1, Y = y / ldDimSamp * diffY + y1;
			LD z = _currentFunc2D(X, Y);
			printf("%.10Lf %.10Lf %.10Lf\n", X, Y, z);
#else
			LD z = _currentFunc2D(x / ldDimSamp * diffX + x1, y / ldDimSamp * diffY + y1);
#endif
			sampling++;
			sum += z;
			var += z * z;
		}
	}

	ssum = var;
	mean = sum / samples;

	error = (x2 - x1) * sqrt(((var / samples) - mean * mean) / samples);

	if(samples < 2)
		variance = 0.0;
	else
		variance = var / (samples - 1) - mean * mean;

	mean *= diffX * diffY;
	return error;
}

inline LD wa(LD x1, LD x2, LD coef1, LD coef2)
{
	LD c1 = coef1 / (coef1 + coef2);
	LD c2 = 1.0 - c1;
	return x1 * c1 + x2 * c2;
}

inline LD was(LD x1, LD x2, LD coef1, LD coef2)
{
	coef1 *= coef1;
	coef2 *= coef2;
	LD c1 = coef1 / (coef1 + coef2);
	LD c2 = 1.0 - c1;
	return x1 * c1 + x2 * c2;
}

/* Adaptivna integracija za jednu dimenziju */
LD adaptiveIntegrate(LD x1, LD x2, LD seg, ULL samples, LD &mean, LD &ssum, LD &variance, int step)
{
	long fsamp = samples / 2;
	long qsamp = fsamp / 2;
	long samp1, samp2;
	LD xm = (x2 - x1) / 2.0;
	LD m1, m2, s1, s2, var1, var2, m11, var11, m22, var22, s11, s22, vsum;
	LD diff = x2 - x1;
	LD frag = (diff / seg);
	LD ovar1, ovar2;

	/* Ako je preostalo jako malo uzoraka ili se doslo do poslednjeg koraka,
	 * onda se radi obicna MC integracija i ne ide se dalje. */
	if(step == 0 || samples < 2)
	{
		mcIntegrate(x1, x2, samples, m1, s1, var1);
		mean = m1;
		variance = var1 * frag;
		ssum = s1;
		return 0.0;
	}
	else
	{
		/* Po cetvrtina uzoraka se rasporedjuje na dve polovine domena */
		mcIntegrate(x1, xm, qsamp, m1, s1, var1);
		mcIntegrate(xm, x2, fsamp - qsamp, m2, s2, var2);

		vsum = var1 + var2;
		ovar1 = var1; ovar2 = var2;

		if(vsum == 0.0)
		{
			variance = (var1 + var2) / 2.0;
			mean = (m1 + m2) / 2.0;
			return 0.0;
		}
		var1 /= vsum;
		var2 /= vsum;

		/* Na osnovu varijanse, druge dve cetvrtine uzoraka se proporcionalno
		 * rasporedjuju. */
		samples -= fsamp;
		fsamp = samples / 2;
		samp1 = var1 * fsamp;
		samp2 = samples - samp1;
		var1 = (LD)samp1 / (LD)samples;
		var2 = (LD)samp2 / (LD)samples;

		/* Funkcija se rekurzivno poziva za polovine domena koristeci izracunatu
		 * raspodelu uzoraka. */
		adaptiveIntegrate(x1, xm, seg, samp1, m11, s11, var11, step - 1);
		adaptiveIntegrate(xm, x2, seg, samp2, m22, s22, var22, step - 1);

		mean = wa(m1, m11, qsamp, samp1) + wa(m2, m22, qsamp, samp2);

		var11 = was(ovar1, var11, 0.5, var1);
		var22 = was(ovar2, var22, 0.5, var2);

		variance = was(var11, var22, 1.0 + var1, 1.0 + var2);

		return 0.0;
	}
}

/* Adaptivna integracija za dve dimenzije, isti princip, samo se domen deli na cetvrtine,
 * a broj uzoraka na osmine. */
LD adaptiveIntegrate2D(LD x1, LD x2, LD y1, LD y2, LD segX, LD segY, ULL samples, LD &mean, LD &ssum, LD &variance, int step)
{
	ULL fsamp = samples / 2;
	ULL qsamp = fsamp / 4;
	ULL samp1, samp2, samp3, samp4;
	LD xm = (x2 + x1) / 2.0;
	LD ym = (y2 + y1) / 2.0;
	LD m1, m2, m3, m4, s1, s2, s3, s4, var1, var2, var3, var4, m11, var11, m22, var22, m33, var33, m44, var44, s11, s22, s33, s44, vsum;
	LD diffX = x2 - x1;
	LD diffY = y2 - y1;
	LD fragX = (diffX / segX);
	LD fragY = (diffY / segY);
	LD ovar1, ovar2, ovar3, ovar4;

	if(step == 0 || samples < 2)
	{
		mcIntegrate2D(x1, x2, y1, y2, samples, m1, s1, var1);
		mean = m1;
		variance = var1 * fragX * fragY;
		ssum = s1;
		return 0.0;
	}
	else
	{
		mcIntegrate2D(x1, xm, y1, ym, qsamp, m1, s1, var1);
		mcIntegrate2D(xm, x2, y1, ym, qsamp, m2, s2, var2);
		mcIntegrate2D(x1, xm, ym, y2, qsamp, m3, s3, var3);
		mcIntegrate2D(xm, x2, ym, y2, qsamp, m4, s4, var4);

		vsum = var1 + var2 + var3 + var4;
		ovar1 = var1; ovar2 = var2;
		ovar3 = var3; ovar4 = var4;

		if(vsum == 0.0)
		{
			variance = (var1 + var2 + var3 + var4) / 4.0;
			mean = (m1 + m2 + m3 + m4) / 4.0;
			return 0.0;
		}
		var1 /= vsum;
		var2 /= vsum;
		var3 /= vsum;
		var4 /= vsum;

		samples -= qsamp * 4;
		fsamp = samples;
		samp1 = var1 * fsamp;
		samp2 = var2 * fsamp;
		samp3 = var3 * fsamp;
		samp4 = fsamp - samp1 - samp2 - samp3;

		var1 = (LD)samp1 / (LD)samples;
		var2 = (LD)samp2 / (LD)samples;
		var3 = (LD)samp3 / (LD)samples;
		var4 = (LD)samp4 / (LD)samples;

		adaptiveIntegrate2D(x1, xm, y1, ym, segX / 2.0, segY / 2.0, samp1, m11, s11, var11, step - 1);
		adaptiveIntegrate2D(xm, x2, y1, ym, segX / 2.0, segY / 2.0, samp2, m22, s22, var22, step - 1);
		adaptiveIntegrate2D(x1, xm, ym, y2, segX / 2.0, segY / 2.0, samp3, m33, s33, var33, step - 1);
		adaptiveIntegrate2D(xm, x2, ym, y2, segX / 2.0, segY / 2.0, samp4, m44, s44, var44, step - 1);

		mean = wa(m1, m11, fsamp, samp1) + wa(m2, m22, fsamp, samp2) + wa(m3, m33, fsamp, samp3) + wa(m4, m44, fsamp, samp4);

		var11 = was(ovar1, var11, 0.5, var1);
		var22 = was(ovar2, var22, 0.5, var2);
		var33 = was(ovar3, var33, 0.5, var3);
		var44 = was(ovar4, var44, 0.5, var4);

		variance = was(was(var11, var22, 1.0 + var1, 1.0 + var2), was(var33, var44, 1.0 + var3, 1.0 + var4), var1 + var2, var3 + var4);

		return 0.0;
	}
}

/* Funkcije omotaci, tu su samo radi jednostavnijeg pristupa spolja */
LD integrateFunction(LD (*mathFunction)(LD), LD lowerLimit, LD upperLimit, IntegrationMethod method,
		ULL samples, LD &mean, LD &variance)
{
	LD temp;
	sampling = 0;
	_currentFunc = mathFunction;

	switch(method)
	{
	case IM_Uniform:
		return uniformIntegrate(lowerLimit, upperLimit, samples, mean, temp, variance);
		break;
	case IM_MonteCarlo:
		return mcIntegrate(lowerLimit, upperLimit, samples, mean, temp, variance);
		break;
	case IM_AdaptiveMC:
		return adaptiveIntegrate(lowerLimit, upperLimit, upperLimit - lowerLimit, samples, mean, temp, variance, _miserSteps);
		break;
	}

	return mean;
}

LD integrateFunction2D(LD (*mathFunction)(LD, LD), LD lowerLimitX, LD upperLimitX, LD lowerLimitY, LD upperLimitY,
		IntegrationMethod method, ULL samples, LD &mean, LD &variance)
{
	LD temp;
	sampling = 0;
	_currentFunc2D = mathFunction;

	switch(method)
	{
	case IM_Uniform:
		return uniformIntegrate2D(lowerLimitX, upperLimitX, lowerLimitY, upperLimitY, samples, mean, temp, variance);
		break;
	case IM_MonteCarlo:
		return mcIntegrate2D(lowerLimitX, upperLimitX, lowerLimitY, upperLimitY, samples, mean, temp, variance);
		break;
	case IM_AdaptiveMC:
		return adaptiveIntegrate2D(lowerLimitX, upperLimitX, lowerLimitY, upperLimitY, upperLimitX - lowerLimitX, upperLimitY - lowerLimitY, samples, mean, temp, variance, _miserSteps);
		break;
	}

	return mean;
}

LD integrateFunction(MAMC_TestCase &test, IntegrationMethod method, ULL samples, LD &mean, LD &variance)
{
	LD err;
	if(test.mathFunction1D != NULL)
	{
		err = integrateFunction(test.mathFunction1D, test.lowerLimitX, test.upperLimitX, method, samples, mean, variance);
		mean *= test.multiplier;
		err *= test.multiplier;
		return err;
	}
	else if(test.mathFunction2D != NULL)
	{
		err = integrateFunction2D(test.mathFunction2D, test.lowerLimitX, test.upperLimitX, test.lowerLimitY, test.upperLimitY, method, samples, mean, variance);
		mean *= test.multiplier;
		err *= test.multiplier;
		return err;
	}
	return 0.0;
}

void setAdaptiveSteps(int steps)
{
	if(steps > 0)
		_miserSteps = steps;
}

ULL getSampleCount()
{
	return sampling;
}

/* Testna funkcije u 1D i 2D varijanti koja ce istu matematicku funkciju da integrise
 * na sva tri podrzana nacina i ispise rezultate. */
void testIntegrations(LD (*func)(LD), LD lowerLimit, LD upperLimit, ULL samples, LD result)
{
	LD m, v, err;
	srand(time(0));

	cout << endl << "-- UNIFORM SAMPLING ---" << endl;
	err = integrateFunction(func, lowerLimit, upperLimit, IM_Uniform, samples, m, v);

	printf("mean : %.10Lf\n", m);
	printf("variance : %.10Lf\n", v);
	printf("delta : %.10Lf\n", abs(m - result));
	printf("est. error : %.10Lf\n", err);
	cout << "samples taken : " << getSampleCount() << endl;

	cout << endl << "-- PLAIN MONTE CARLO ---" << endl;
	err = integrateFunction(func, lowerLimit, upperLimit, IM_MonteCarlo, samples, m, v);

	printf("mean : %.10Lf\n", m);
	printf("variance : %.10Lf\n", v);
	printf("delta : %.10Lf\n", abs(m - result));
	printf("est. error : %.10Lf\n", err);
	cout << "samples taken : " << getSampleCount() << endl;


	cout << endl << "-- ADAPTIVE MONTE CARLO ---" << endl;
	err = integrateFunction(func, lowerLimit, upperLimit, IM_AdaptiveMC, samples, m, v);

	printf("mean : %.10Lf\n", m);
	printf("variance : %.10Lf\n", v);
	printf("delta : %.10Lf\n", abs(m - result));
	printf("est. error : %.10Lf\n", err);
	cout << "samples taken : " << getSampleCount() << endl;
}

void testIntegrations2D(LD (*func)(LD, LD), LD lowerLimitX, LD upperLimitX, LD lowerLimitY, LD upperLimitY, ULL samples, LD result)
{
	LD m, v, err;
	srand(time(0));

	cout << endl << "-- UNIFORM SAMPLING ---" << endl;
	err = integrateFunction2D(func, lowerLimitX, upperLimitX, lowerLimitY, upperLimitY, IM_Uniform, samples, m, v);

	printf("mean : %.10Lf\n", m);
	printf("variance : %.10Lf\n", v);
	printf("delta : %.10Lf\n", abs(m - result));
	printf("est. error : %.10Lf\n", err);
	cout << "samples taken : " << getSampleCount() << endl;

	cout << endl << "-- PLAIN MONTE CARLO ---" << endl;
	err = integrateFunction2D(func, lowerLimitX, upperLimitX, lowerLimitY, upperLimitY, IM_MonteCarlo, samples, m, v);

	printf("mean : %.10Lf\n", m);
	printf("variance : %.10Lf\n", v);
	printf("delta : %.10Lf\n", abs(m - result));
	printf("est. error : %.10Lf\n", err);
	cout << "samples taken : " << getSampleCount() << endl;


	cout << endl << "-- ADAPTIVE MONTE CARLO ---" << endl;
	err = integrateFunction2D(func, lowerLimitX, upperLimitX, lowerLimitY, upperLimitY, IM_AdaptiveMC, samples, m, v);

	printf("mean : %.10Lf\n", m);
	printf("variance : %.10Lf\n", v);
	printf("delta : %.10Lf\n", abs(m - result));
	printf("est. error : %.10Lf\n", err);
	cout << "samples taken : " << getSampleCount() << endl;
}

void testIntegrations(MAMC_TestCase &test, ULL samples)
{
	if(test.mathFunction1D != NULL)
		testIntegrations(test.mathFunction1D, test.lowerLimitX, test.upperLimitX, samples, test.result);
	else if(test.mathFunction2D != NULL)
		testIntegrations2D(test.mathFunction2D, test.lowerLimitX, test.upperLimitX, test.lowerLimitY, test.upperLimitY, samples, test.result);
}
